/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef itkBioCellularAggregate_hxx
#define itkBioCellularAggregate_hxx

#include "itkBioCellularAggregate.h"

namespace itk
{
namespace bio
{
template <unsigned int NSpaceDimension>
CellularAggregate<NSpaceDimension>::CellularAggregate()
{
  typename BioCellType::ColorType color;
  color.SetRed(1.0);
  color.SetGreen(1.0);
  color.SetBlue(1.0);
  BioCellType::SetDefaultColor(color);

  m_Mesh = MeshType::New();

  m_Mesh->SetCellsAllocationMethod(MeshEnums::MeshClassCellsAllocationMethod::CellsAllocatedDynamicallyCellByCell);
  m_Mesh->SetPoints(PointsContainer::New());
  m_Mesh->SetPointData(PointDataContainer::New());
  m_Mesh->SetCells(VoronoiRegionsContainer::New());

  m_Iteration = 0;
  m_ClosestPointComputationInterval = 5;

  m_FrictionForce = 1.0f;
}

template <unsigned int NSpaceDimension>
CellularAggregate<NSpaceDimension>::~CellularAggregate()
{
  this->KillAll();
}

template <unsigned int NSpaceDimension>
unsigned int
CellularAggregate<NSpaceDimension>::GetNumberOfCells() const
{
  return m_Mesh->GetPointData()->Size();
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::SetGrowthRadiusLimit(double value)
{
  BioCellType::SetGrowthRadiusLimit(value);
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::SetGrowthRadiusIncrement(double value)
{
  BioCellType::SetGrowthRadiusIncrement(value);
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::PrintSelf(std::ostream & os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << "Cellular aggregate " << m_Mesh << std::endl;
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::Remove(CellBase * cellbase)
{
  auto * cell = dynamic_cast<BioCellType *>(cellbase);

  if (!cell)
  {
    itk::ExceptionObject exception;
    exception.SetDescription("Trying to remove a null pointer to cell");
    exception.SetLocation("CellularAggregate::Remove(BioCellType*)");
    throw exception;
  }

  IdentifierType id = cell->GetSelfIdentifier();

  typename MeshType::CellAutoPointer region;
  bool                               regionExist = m_Mesh->GetCell(id, region);
  if (regionExist)
  {
    auto * realRegion = dynamic_cast<VoronoiRegionType *>(region.GetPointer());

    if (!realRegion)
    {
      itkExceptionMacro("CellularAggregate::Remove() couldn't dynamic_cast region ");
    }
    else
    {
      //
      // Notify all the neighbors that this cell is going away
      //
      typename VoronoiRegionType::PointIdIterator neighbor = realRegion->PointIdsBegin();
      typename VoronoiRegionType::PointIdIterator end = realRegion->PointIdsEnd();
      while (neighbor != end)
      {
        const IdentifierType               neighborId = (*neighbor);
        typename MeshType::CellAutoPointer cellPointer;
        bool                               neighborVoronoiExist = m_Mesh->GetCell(neighborId, cellPointer);
        if (neighborVoronoiExist)
        {
          auto * vregion = dynamic_cast<VoronoiRegionType *>(cellPointer.GetPointer());
          if (!vregion)
          {
            std::cerr << "CellularAggregate::Add() Failed to find a region" << std::endl;
          }
          else
          {
            vregion->RemovePointId(id);
          }
        }
        neighbor++;
      }

      // We can now remove the entry from the list of cells in the Mesh
      m_Mesh->GetCells()->DeleteIndex(id);

      // and then release the memory used by the VoronoiRegion
      delete realRegion;
    }
  }
  else
  {
    itkExceptionMacro(" Region " << id << " doesn't exist ");
  }

  m_Mesh->GetPoints()->DeleteIndex(id);
  m_Mesh->GetPointData()->DeleteIndex(id);

  // Finally we can delete the BioCell;
  delete cell;
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::GetVoronoi(IdentifierType cellId, VoronoiRegionAutoPointer & voronoiPointer) const
{
  typename MeshType::CellAutoPointer cellPointer;

  bool voronoiExists = m_Mesh->GetCell(cellId, cellPointer);
  if (!voronoiExists)
  {
    itk::ExceptionObject exception;
    exception.SetDescription("voronoi region does not exist in the container");
    exception.SetLocation("GetVoronoi( IdentifierType");
    throw exception;
  }

  auto * region = dynamic_cast<VoronoiRegionType *>(cellPointer.GetPointer());

  voronoiPointer.TakeNoOwnership(region);
  if (cellPointer.IsOwner())
  {
    voronoiPointer.TakeOwnership();
    cellPointer.ReleaseOwnership();
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::SetEgg(BioCellType * cell, const PointType & position)
{
  VectorType perturbation = position.GetVectorFromOrigin();

  this->Add(cell, perturbation);
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::Add(CellBase * cell)
{
  VectorType perturbation;

  perturbation.Fill(0.0);
  this->Add(cell, perturbation);
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::Add(CellBase * cellA, CellBase * cellB, double perturbationLength)
{
  // Create a perturbation for separating the daugther cells
  typename BioCellType::VectorType perturbationVector;
  for (unsigned int d = 0; d < NSpaceDimension; d++)
  {
    perturbationVector[d] = vnl_sample_uniform(-1.0f, 1.0f);
  }

  const double norm = perturbationVector.GetNorm();
  if (itk::Math::abs(norm) > 1e-10)
  {
    perturbationVector *= perturbationLength / norm;
  }
  else
  {
    // this event should rarely happen... very rarely
    perturbationVector[0] = perturbationLength;
  }

  this->Add(cellA, perturbationVector);
  this->Add(cellB, -perturbationVector);

  const IdentifierType cellAId = cellA->GetSelfIdentifier();
  const IdentifierType cellBId = cellB->GetSelfIdentifier();

  VoronoiRegionAutoPointer cellAptr;
  VoronoiRegionAutoPointer cellBptr;

  this->GetVoronoi(cellAId, cellAptr);
  this->GetVoronoi(cellBId, cellBptr);

  cellAptr->AddPointId(cellBId);
  cellBptr->AddPointId(cellAId);
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::Add(CellBase * cellBase, const VectorType & perturbation)
{
  auto * cell = dynamic_cast<BioCellType *>(cellBase);
  if (cell == nullptr)
  {
    itkExceptionMacro(<< "dynamic_cast failed.");
  }
  IdentifierType newcellId = cell->GetSelfIdentifier();
  IdentifierType newcellparentId = cell->GetParentIdentifier();

  PointType       position;
  CellAutoPointer selfVoronoi;

  // If the cell does not have a parent
  // from which receive a position
  if (!newcellparentId)
  {
    position.Fill(0.0);
    selfVoronoi.TakeOwnership(new VoronoiRegionType);
  }
  else
  {
    bool parentPositionExist = m_Mesh->GetPoint(newcellparentId, &position);
    if (!parentPositionExist)
    {
      itk::ExceptionObject exception;
      exception.SetDescription("Parent cell does not exist in the container");
      exception.SetLocation("CellularAggregate::Add( cell * )");
      throw exception;
    }

    VoronoiRegionAutoPointer parentVoronoi;
    this->GetVoronoi(newcellparentId, parentVoronoi);
    parentVoronoi->MakeCopy(selfVoronoi);
  }

  position += perturbation;

  m_Mesh->SetCell(newcellId, selfVoronoi);
  m_Mesh->SetPoint(newcellId, position);
  m_Mesh->SetPointData(newcellId, cell);

  cell->SetCellularAggregate(this);

  // Add this new cell as neighbor to cells in its neighborhood

  typename VoronoiRegionType::PointIdIterator neighbor = selfVoronoi->PointIdsBegin();
  typename VoronoiRegionType::PointIdIterator end = selfVoronoi->PointIdsEnd();

  while (neighbor != end)
  {
    const IdentifierType               neighborId = (*neighbor);
    typename MeshType::CellAutoPointer cellPointer;
    bool                               neighborVoronoiExist = m_Mesh->GetCell(neighborId, cellPointer);
    if (neighborVoronoiExist)
    {
      auto * region = dynamic_cast<VoronoiRegionType *>(cellPointer.GetPointer());

      if (!region)
      {
        std::cerr << "CellularAggregate::Add() Failed to find a region" << std::endl;
      }
      else
      {
        region->AddPointId(newcellId);
      }
    }
    neighbor++;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::AdvanceTimeStep()
{
  if (m_Iteration % m_ClosestPointComputationInterval == 0)
  {
    this->ComputeClosestPoints();
  }

  this->ComputeForces();
  this->UpdatePositions();

  CellsIterator begin = m_Mesh->GetPointData()->Begin();
  CellsIterator end = m_Mesh->GetPointData()->End();

  CellsIterator cell = begin;

  while (cell != end)
  {
    BioCellType * theCell = cell.Value();
    theCell->AdvanceTimeStep();
    ++cell;
    if (theCell->MarkedForRemoval())
    {
      this->Remove(theCell);
    }
  }

  this->InvokeEvent(IterationEvent());

  m_Iteration++;
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::KillAll()
{
  if (!m_Mesh)
  {
    return;
  }

  CellsIterator cell = m_Mesh->GetPointData()->Begin();
  CellsIterator end = m_Mesh->GetPointData()->End();

  while (cell != end)
  {
    delete (cell.Value());
    ++cell;
  }

  BioCellType::ResetCounter();
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::ClearForces()
{
  CellsIterator cell = m_Mesh->GetPointData()->Begin();
  CellsIterator end = m_Mesh->GetPointData()->End();

  while (cell != end)
  {
    cell.Value()->ClearForce();
    ++cell;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::UpdatePositions()
{
  CellsConstIterator cellIt = m_Mesh->GetPointData()->Begin();
  CellsConstIterator end = m_Mesh->GetPointData()->End();

  PointType position;

  position.Fill(0);

  while (cellIt != end)
  {
    BioCellType *  cell = cellIt.Value();
    IdentifierType cellId = cell->GetSelfIdentifier();
    m_Mesh->GetPoint(cellId, &position);
    const VectorType force = cell->GetForce();
    if (force.GetNorm() > m_FrictionForce)
    {
      position += force / 50.0;
    }
    m_Mesh->SetPoint(cellId, position);
    cellIt++;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::ComputeForces()
{
  // Clear all the force accumulators
  this->ClearForces();

  CellsConstIterator cell1It = m_Mesh->GetPointData()->Begin();
  CellsConstIterator end = m_Mesh->GetPointData()->End();

  // compute forces
  while (cell1It != end)
  {
    const IdentifierType cell1Id = cell1It.Index();

    BioCellType * cell1 = cell1It.Value();

    PointType position1;
    position1.Fill(0);
    m_Mesh->GetPoint(cell1Id, &position1);

    const double rA = cell1->GetRadius();

    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1Id, voronoiRegion);

    typename VoronoiRegionType::PointIdIterator neighbor = voronoiRegion->PointIdsBegin();
    typename VoronoiRegionType::PointIdIterator vend = voronoiRegion->PointIdsEnd();

    while (neighbor != vend)
    {
      const IdentifierType cell2Id = (*neighbor);

      BioCellType * cell2 = nullptr;
      PointType     position2;

      if (!m_Mesh->GetPoint(cell2Id, &position2))
      {
        neighbor++; // if the neigbor has been removed, skip it
        continue;
      }
      m_Mesh->GetPointData(cell2Id, &cell2);

      const double rB = cell2->GetRadius();

      typename BioCellType::VectorType relativePosition = position1 - position2;

      const double distance = relativePosition.GetNorm();

      if (distance < (rA + rB) / 2.0)
      {
        const double                     factor = 2.0 * BioCellType::GetGrowthRadiusLimit() / distance;
        typename BioCellType::VectorType force = relativePosition * factor;
        cell1->AddForce(force);
        cell2->AddForce(-force);
      }
      else if (distance < rA + rB)
      {
        typename BioCellType::VectorType force = relativePosition;
        cell1->AddForce(force);
        cell2->AddForce(-force);
      }

      neighbor++;
    }

    cell1It++;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::ComputeClosestPoints()
{
  PointsConstIterator beginPoints = m_Mesh->GetPoints()->Begin();
  PointsConstIterator endPoints = m_Mesh->GetPoints()->End();

  PointsConstIterator point1It = beginPoints;

  while (point1It != endPoints)
  {
    PointType position1 = point1It.Value();

    PointsConstIterator point2It = beginPoints;

    BioCellType * cell1 = nullptr;

    IdentifierType cell1Id = point1It.Index();
    m_Mesh->GetPointData(cell1Id, &cell1);
    const double radius = cell1->GetRadius();
    const double limitDistance = radius * 4.0;

    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1Id, voronoiRegion);

    voronoiRegion->ClearPoints();

    while (point2It != endPoints)
    {
      if (point2It == point1It)
      {
        point2It++;
        continue;
      }

      const PointType & position2 = point2It.Value();

      typename BioCellType::VectorType relativePosition = position1 - position2;

      const double distance = relativePosition.GetNorm();
      if (distance < limitDistance)
      {
        voronoiRegion->AddPointId(point2It.Index());
      }
      point2It++;
    }

    point1It++;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::DumpContent(std::ostream & os) const
{
  CellsConstIterator beginCell = m_Mesh->GetPointData()->Begin();
  CellsConstIterator endCell = m_Mesh->GetPointData()->End();

  CellsConstIterator cell1It = beginCell;

  os << "Cell Identifiers " << std::endl;
  while (cell1It != endCell)
  {
    os << cell1It.Index() << " == ";
    os << cell1It.Value()->GetSelfIdentifier() << std::endl;
    cell1It++;
  }

  os << std::endl << "Points " << std::endl;
  PointsConstIterator pointIt = m_Mesh->GetPoints()->Begin();
  PointsConstIterator endPoint = m_Mesh->GetPoints()->End();

  while (pointIt != endPoint)
  {
    os << pointIt.Index() << "  ";
    os << pointIt.Value() << std::endl;
    pointIt++;
  }

  os << std::endl << "Neighborhoods " << std::endl;
  cell1It = beginCell;
  while (cell1It != endCell)
  {
    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1It.Index(), voronoiRegion);
    typename VoronoiRegionType::PointIdIterator neighbor = voronoiRegion->PointIdsBegin();
    typename VoronoiRegionType::PointIdIterator end = voronoiRegion->PointIdsEnd();

    while (neighbor != end)
    {
      os << (*neighbor) << "   ";
      neighbor++;
    }
    os << std::endl;

    cell1It++;
  }
}

template <unsigned int NSpaceDimension>
void
CellularAggregate<NSpaceDimension>::AddSubstrate(SubstrateType * substrate)
{
  SubstratePointer smartPointer(substrate);

  m_Substrates.push_back(smartPointer);
}

template <unsigned int NSpaceDimension>
typename CellularAggregate<NSpaceDimension>::SubstratesVector &
CellularAggregate<NSpaceDimension>::GetSubstrates()
{
  return m_Substrates;
}

template <unsigned int NSpaceDimension>
typename CellularAggregate<NSpaceDimension>::SubstrateValueType
CellularAggregate<NSpaceDimension>::GetSubstrateValue(IdentifierType cellId, unsigned int substrateId) const
{
  PointType cellPosition;
  bool      cellPositionExists = m_Mesh->GetPoint(cellId, &cellPosition);

  if (!cellPositionExists)
  {
    std::cerr << " Cell position doesn't exist for cell Id = ";
    std::cerr << cellId << std::endl;
    return itk::NumericTraits<SubstrateValueType>::ZeroValue();
  }

  SubstratePointer substrate = m_Substrates[substrateId];

  typename SubstrateType::IndexType index;

  substrate->TransformPhysicalPointToIndex(cellPosition, index);

  SubstrateValueType value = 0;
  if (substrate->GetBufferedRegion().IsInside(index))
  {
    value = substrate->GetPixel(index);
  }

  return value;
}
} // end namespace bio
} // end namespace itk

#endif
